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Abstract. For axially symmetric even-even nuclei, we solve the Hartree-Fock-Bogoliubov 
(HFB) equations on a 2-D grid in cylindrical coordinates. The Skyrme SLy4 interaction is 
used for the mean field and a zero-range interaction for the pairing field. After decoupling the 
HFB lattice equations, we obtain quasiparticle states with equivalent single-particle energies 
up to f00 MeV or more. We present results for the Krypton isotope chain up to the two- 
^ ' neutron dripline, including two-neutron separation energies, pairing gaps, and quadrupole 

, deformations for the ground states and isomeric minima. 
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g ; 1 INTRODUCTION 

Experimental studies of the limits of nuclear stability provide stringent tests 
for nuclear structure theories. While nuclei near the proton dripline have 
been well explored, little is known about the limits of nuclear binding on 
the neutron-rich side and about the exact location of the neutron dripline. 
Only for the eight lightest nuclei has the neutron dripline been reached, and 
the nuclear chart exhibits thousands of nuclear isotopes still to be examined 
with the next generation of Radioactive Ion Beam facilities. Another limit 
to nuclear stability is the superheavy element region which is formed by a 
delicate balance between strong Coulomb repulsion and additional binding 
due to closed shells. 

For very light nuclei with A < 12, ab initio calculations have been car- 
ried out using nucleon-nucleon and three- nucleon interactions. For heavier 
nuclei, approximations are required; the most widely used approach is the 
self-consistent mean field method. Both non-relativistic mean field theories 
[H El El H El El E] and relativistic versions [El El EI] have been developed. 

In the vicinity of the driplines, pairing correlations increase dramati- 
cally and it is essential to treat both the mean field and the pairing field 
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self-consistently within the Hartree-Fock-Bogoliubov (HFB) formalism pQ. In 
this region, not only does one have to consider "well-bound" single-particle 
states but also "weakly-bound" states with large spatial extent, giving rise to 
neutron skins. All of these features represent major challenges for the numer- 
ical solution. While the HF(B) theories describe the ground state properties 
of nuclei, their excited states can be obtained with the (quasiparticle) random 
phase approximation (Q)RPA [TT1 IT21 IT5] . 

Traditionally, the HFB equations have been solved by expanding the 
quasiparticle wavefunctions in a truncated harmonic oscillator basis [2] . This 
works very well near the line of /^-stability because only "well-bound" states 
need to be considered. However, as one approaches the driplines, the numerical 
solution becomes more challenging: in practice, it is very difficult to represent 
continuum states as superpositions of bound harmonic oscillator states because 
the former show oscillatory behavior at large distances while the latter decay 
exponentially. In this case, coordinate-space representations have a distinct 
advantage because "bound" and "continuum" states (in the box) are treated 
on an equal footing, with no bias towards bound states. 



2 HFB EQUATIONS IN COORDINATE SPACE 



Over a period of several years, our research group has developed the HFB- 
2D-LATTICE code, written in Fortran 95, which solves the HFB equations 
for deformed, axially symmetric even-even nuclei in coordinate space on a 2- 
D lattice [El EJ QH QH QU [20] . The novel feature of our HFB code is that 
it generates high-energy continuum states with an equivalent single-particle 
energy of more than 100 MeV. QRPA calculations require the inclusion of 
continuum states with such high energy for the description of collective giant 
multipole resonances. 

A detailed description of our theoretical method has been published in 
Ref. [6]; in the following, we give a brief summary. In coordinate space repre- 
sentation, the HFB Hamiltonian and the quasiparticle wavefunctions depend 
on the distance vector r, spin projection a = ±|, and isospin projection 
q = ±^ (corresponding to protons and neutrons, respectively). In the HFB 
formalism, there are two types of quasiparticle wavefunctions, <f>i and fa, which 
are bi-spinors of the form 



Their dependence on the quasiparticle energy is denoted by the index a 
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for simplicity. The quasiparticle wavefunctions determine the normal density 



A?0) 

and the pairing density 

Pg( r ) = 
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El>0a=- 
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In the present work, we use the Skyrme SLy4 effective N-N interaction [2T] 
for the mean field (p-h and h-p channels), and a zero-range pairing force (p-p 
and h-h channels). The pairing strength has been adjusted to reproduce the 
measured average neutron paring gap of 1.245 MeV in 120 Sn. For these types 
of effective interactions, the particle mean field Hamiltonian h and the pairing 
field Hamiltonian h are diagonal in isospin space and local in position space 



h(raq,r'a'q') = 8 q>q < 5{r - r')h q ,(r) 



and 



h(raq, r a q 



J q,q' 
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In spin-space, the mean field Hamiltonian is represented by the 2x2 matrix 



h q (r) 




h q n (r) 
h q u (r) 



(6) 



which operates on the bi-spinor wavefunctions in eq. The pairing field 
Hamiltonian h has a similar mathematical structure. The HFB equations 
have the following structure in spin-space [6]: 





E q 

J —'r\i 



(7) 



where X q denotes the Fermi level. From now on, we drop the isospin label q, 
for simplicity. The HFB equations may be recast in the form 



Hcf> a 

with the four-spinor wavefunctions 



(8) 
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and the quasiparticle "super-Hamiltonian" (4x4 matrix in spin space) is given 

This quasiparticle Hamiltonian has the well-known property [221 23] that for 
every eigenvector </>+ = (U a ,V a ) with positive eigenvalue E+, there is a cor- 
responding eigenvector <f>~ = (V*,U*) with negative eigenvalue E~ = —E+, 
respectively 

H<£ = E±<i>±. (11) 

It is forbidden to choose the E+ and E~ solutions at the same time, otherwise 
it is impossible to fulfill the anti-commutation relations for the quasiparticle 
operators [22]. The quasiparticle energy spectrum is discrete for \E a \ < — A 
and continuous for \E a \ > —A. For even-even nuclei it is customary to solve 
the HFB equations for positive quasiparticle energies and consider all negative 
energy states as occupied in the HFB ground state pQ. 



3 2-D LATTICE REPRESENTATION 

In the case of even-even nuclei with axially symmetric shapes, we represent 
the HFB equations in cylindrical coordinates ((f), r, z) and eliminate the depen- 
dence on the coordinate (j>, resulting in a reduced 2-D problem. For a given 
angular momentum projection quantum number 0, we solve the HFB equa- 
tions on a 2-D grid {r a ,za) where a = l,...,m and (3 = l,...,n. In practice, 
angular momentum projections $7 = 1/2, 3/2, 21/2 are taken into account. 
Typically, in radial (r) direction, the lattice extends from — 15 fm, and 
in symmetry axis (z) direction from — 15, +15 fm, with a lattice spacing 
of about 0.8 fm in the central region. For the lattice representation of the 
HFB Hamiltonian we use a hybrid method in which derivative operators are 
constructed using the Galerkin method [E]; this amounts to a global error 
reduction. Local potentials are represented by the Basis-Spline collocation 
method [25j[26l[27] (local error reduction). 

On the lattice, the four-spinor wavefunction in coordinate space 4>^ (r, z) 
becomes an array 4>(N) of length N = 4 x m x n and the HFB Hamiltonian 
is transformed into a N x N matrix: 

Y j H;$ = E^ Ou = l,...,iV) . (12) 

v=l 

The diagonalization of the HFB Hamiltonian using LAPACK routines yields 
quasiparticle states with energies up to 100 MeV or more. In the calculations 
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of observables, all quasiparticles states with equivalent single particle energy 
up to the cut-off e max = 60 MeV are taken into account. 

Production runs of our HFB code are carried out on an IBM-SP massively 
parallel supercomputer and on a local LINUX cluster using OPENMP/MPI 
message passing. Parallelization is possible for different angular momentum 
states Q and isospins (p/n). 



4 DECOUPLED HFB EQUATIONS 

The HFB equations ([7]) are coupled, i.e. the HFB Hamiltonian in eq. (I10p 
mixes the quasiparticle states <^>i jQ and (f>2,a- We generalize the recipe given in 
Ref. [24J to decouple the HFB equations which results in a smaller diagonal- 
ization problem. 

This is accomplished by the unitary transformation 

U - BZ - * = 7f(!-i)' z = 7f(-i0' (13) 

First we apply the complex unitary transformation Z to the quasiparticle 
Hamiltonian resulting in a transformed Hamiltonian in which the pairing fields 
on the diagonal appear with opposite sign 

* = ^' = (i, + -f )• <") 

where we have introduced the abbreviation 

h' = h-X. (15) 

Next we apply the real transformation R and obtain a new Hamiltonian Hjj 
in block off-diagonal form 

H U = UHU- 1 =( j ~ - ih ' + ~ h ) (i 6 ) 

^ +ih' + h J v ' 

with the transformed eigenstates 

xi = U4>i (17) 

resulting in 

Huxt = Etxt ■ (18) 

Operating from the left with the transformed Hamiltonian Hu onto eq. (|18|) 
we obtain 

Hfa'c = E 2 J a ■ (19) 
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It turns out that the operator H^j has the desired block-diagonal form 
_ / (-M + h)(+ih' + h) \ 

Hu ~ { o (+ih' + h)(-ih' + h) ) • ( j 

The new eigenstates x a are > * n general, linear combinations of the negative and 
positive states Xa because we lost information about the sign of the energy 
E a by squaring the Hamiltonian 

Xa = W + xi +W~Xa ( 21 ) 

where w + and w~ are complex constants. 

Prom the structure of the Hamiltonian Hy in eq. (|20p we see that the 
eigenstates x'a l an d x'a.% m ec L- p9j) are now decoupled. Therefore, we only 
need to diagonalize the upper and lower blocks of the Hamiltonian Hfj sep- 
arately, resulting in a speed-up by nearly a factor of 4 in the numerical cal- 
culations compared to the original coupled diagonalization problem, eq. fl7|). 
Further details regarding the numerical solution of the decoupled HFB equa- 
tions can be found in Ref. 



5 RESULTS FOR THE KRYPTON ISOTOPES 

Recently we have studied the ground state properties of the neutron-rich Zirco- 
nium isotopes (Z = 40) up to the two-neutron dripline [TBI ED]- In particular, 
we have calculated binding energies, two-neutron separation energies, normal 
densities and pairing densities, mean square radii, quadrupole moments, and 
pairing gaps. In Ref. pjj], two theoretical approaches were presented and 
compared: a) our HFB-2D-LATTICE code, and b) an expansion in a trans- 
formed harmonic oscillator basis with 20 shells (HFB-2D-THO) by Stoitsov et 
al. [HE]. In general, both theoretical methods were found to be in excellent 
agreement for the entire isotope chain. In particular, 122 Zr is predicted to be 
the dripline nucleus. 

In this section, we present our HFB-2D-LATTICE code results for the 
ground state properties of the neutron-rich even-even Krypton isotopes {Z = 
36) in the mass region A = 104 — 118. Fig. [ljshows the two-neutron separation 
energies obtained from the calculated binding energies. Experimental data for 
the Krypton isotopes calculated in this work are not available because of the 
large N/Z ratios > 1.8. Compared to the two-neutron separation energies for 
the Zr isotope chain for the same mass number (see Fig. 1 of Ref. [19J) the 
separation energy values are only about half the size. Our HFB calculations 
predict that 116 Kr is the last stable nucleus against two-neutron emission. 

Like the Zirconium isotopes studied in Ref. [19J, the Krypton isotope 
chain reveals a rich variety of shapes. In particular, we observe a transition 



6 




from prolate to oblate shapes and back to prolate shapes in the mass re- 
gion A = 104 — 118. Specifically, the quadrupole deformations decrease from 
/?2 =0.38 ( 104 Kr) through /?2=-0.31 ( 108 Kr) and then increase again until we 
reach (32=0.1 for 114 ' 116 Kr. A similar trend for the shape transitions has been 
found by two other models, the Finite Range Droplet Model calculations |28j 
and by the Relativistic Mean-Field calculations [9]. These latter two models 
use a simple BCS pairing interaction which is not reliable for nuclei in the 
vicinity of the two-neutron dripline. 

Table 1: Shape coexistence in the Kr isotopes. The first column lists properties 
of the ground state minimum (quadrupole deformation and binding energy), 
and subsequent columns give results for shape isomeric states . 



isotope 


1st min. 


2nd min. 


3rd min. 






E[MeV] 


/3 2 n 




E[MeV] 


0$ 




E[MeV] 


lU4 Kr 


0.38 


0.35 


-821.41 


0.09 


0.1 


-816.44 


-0.26 


-0.25 


-821.30 


lU6 Kr 


0.36 


0.35 


-825.65 








-0.29 


-0.28 


-825.29 


liu Kr 


-0.15 


-0.15 


-832.14 


0.07 


0.08 


-830.45 


0.40 


0.36 


-829.94 
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We also find evidence for shape coexistence in several neutron-rich Kryp- 
ton isotopes. In Table ([1]) we list the quadrupole deformation and binding 
energy of the ground state PES minimum and several shape-isomeric minima 
for the isotopes I04,i06,li0j^ r _ -^y e no ti ce that the energy difference between 
the prolate and oblate minima for 104 > 106 Kr are very small (of order of 10-40 
keV). 

A plot of the rms-radii for the Krypton isotope chain [2U] shows the 
development of a neutron skin as evidenced by the large difference between 
proton and neutron rms-radii. A sudden rms-radius drop observed at A=108- 
110 coincides with the nuclear shape transition from prolate to oblate in this 
region. A neutron skin development has also been seen in our HFB lattice 
calculations of the Sulfur isotopes (Fig. 4 in Ref . [16] ) and in the Zr isotopes 
(Fig. 4 in Ref. US]). 

1| — 1 — I — 1 — I — 1 — I — 1 — I — 1 — I — 1 — I — 1 — I — 1 — I 




6 68 70 72 74 76 78 80 82 

N 

Figure 2: Average neutron pairing gap for the chain of Krypton isotopes. 

In general, the pairing densities and average pairing gaps for protons and 
neutrons are the two observables which are most sensitive to the properties 
of the continuum states. A study of the spectral distribution of the pairing 
density for the Zirconium and Krypton isotopes [17} [20] shows that it is peaked 
at the Fermi energy A and reaches high into the continuum for neutron-rich 
nuclei. In Figures [2] and [3] we show the average pairing gaps for neutrons and 
protons, respectively. 

We observe that the Krypton pairing gaps in the A = 104 — 118 mass 
region vary between 0.1 and 0.9 MeV. Starting at neutron number N = 68, 
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Figure 3: Average proton pairing gap for the chain of Krypton isotopes. 



the neutron pairing gap increases from 0.1 MeV until it reaches a maximum 
value of 0.9 MeV at N = 76; for higher neutron numbers it decreases quickly 
back to 0.1 MeV. By contrast, the proton pairing gap shows an oscillatory 
structure in the same mass region. 
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